# Convenient jupyter setup
%load_ext autoreload
%autoreload 2
#%load_ext rpy2.ipython

import re
import pathlib
import geopandas as gpd
import pandas as pd
import rioxarray as rxr

from rasterio.enums import Resampling

# Visualisation imports
from holoviews import opts
import holoviews as hv
from holoviews.operation.datashader import rasterize
from holoviews.element import tiles
hv.extension("bokeh")

import numpy as np
from src.constants import (PAISAGENSLIDAR_PATH, EBALIDAR_PATH, DATA_PATH, 
                           WGS84, WEBMERCATOR, SIRGAS_BRAZIL)

from src.utils.os import list_content, latest_matching_file
from src.visualisation.image_loader import fetch_metrics, visualise, load_image, load_raster, load_bounds, VISUALISATION_CACHE, stylise, create_raster_image
from src.visualisation.visualisation_params import *

# Global parameter whether to save
save = False
esri = tiles.ESRI().redim(x="Easting", y="Northing")

from src.data.jrc_loading import get_jrc_data
from src.processing.canopy_height_models import longo2016_acd, asner2014_acd, asner2014_colombia_acd

Load master data

paisagens_master = gpd.read_feather(PAISAGENSLIDAR_PATH / "paisagens_master_table.feather")
eba_master = gpd.read_feather(EBALIDAR_PATH / "eba_master_table.feather")

#paisagens_master.iloc[0]

Define sample

paisagens_survey_name = "SAN_A02_2014"
eba_survey_name = "NP_T-0859"
eba_survey_name2 = "NP_T-0408"

p_sample = paisagens_master.loc[paisagens_survey_name]
e_sample = eba_master.loc[eba_survey_name]
e_sample2 = eba_master.loc[eba_survey_name2]

samples = [p_sample, e_sample, e_sample2]

lon_min = min([sample.lon_min for sample in samples])
lat_min = min([sample.lat_min for sample in samples])
lon_max = max([sample.lon_max for sample in samples])
lat_max = min([sample.lat_max for sample in samples])

Get biomass data

# Get bounds
l,b,r,t = gpd.GeoSeries(p_sample.geometry_bbox, crs=SIRGAS_BRAZIL).to_crs(WEBMERCATOR).bounds.values[0]

#xs, ys = np.array(list(gpd.GeoSeries(sample.geometry_metrics, crs=SIRGAS_BRAZIL).to_crs(WEBMERCATOR)[0].exterior.coords)).T
#poly = hv.Polygons({"x": xs, "y": ys}).redim(x="Easting", y="Northing").opts(alpha=0.2)

resolution = 50
chms = [load_raster(sample.name, "canopy_height_model", 1, from_crs=sample.crs, to_crs=None) for sample in samples]
chm_means = [chm.rio.reproject(chm.rio.crs, resolution=resolution, resampling=Resampling.average) for chm in chms]

# Filter high outliers
chm_means = [chm.where(chm < 100) for chm in chm_means]

longo2016_acds = [longo2016_acd(chm_mean.where(chm_mean < 100)) for chm_mean in chm_means]
asner2014_acds = [asner2014_acd(chm_mean.where(chm_mean < 100)) for chm_mean in chm_means]
asner2014_colombia_acds = [asner2014_colombia_acd(chm_mean.where(chm_mean < 100)) for chm_mean in chm_means]

Load JRC data

transition_map = get_jrc_data(lon_min, lat_min, lon_max, lat_max, dataset="TransitionMap_MainClasses").squeeze().rio.reproject(WEBMERCATOR)
first_deforestation = get_jrc_data(lon_min, lat_min, lon_max, lat_max, dataset="DeforestationYear").squeeze().rio.reproject(WEBMERCATOR)
first_degradation = get_jrc_data(lon_min, lat_min, lon_max, lat_max, dataset="DegradationYear").squeeze().rio.reproject(WEBMERCATOR)
last_deforestation = get_jrc_data(lon_min, lat_min, lon_max, lat_max, dataset="LastDeforestationYear").squeeze().rio.reproject(WEBMERCATOR)
last_degradation = get_jrc_data(lon_min, lat_min, lon_max, lat_max, dataset="LastDegradationYear").squeeze().rio.reproject(WEBMERCATOR)

(visualise("SAN_A02_2014", "TransitionMap", 30, transition_map) 
 + visualise(p_sample.name, "canopy_height_model", 1, datashader_aggregator="mean", from_crs=p_sample.crs)
 + visualise("SAN_A02_2014", "LastDeforestation", 30, last_deforestation.where(last_deforestation > 1) )
 + visualise("SAN_A02_2014", "FirstDeforestation", 30, first_deforestation.where(first_deforestation > 1))
 + visualise("SAN_A02_2014", "LastDegradation", 30, last_degradation.where(last_degradation > 1) )
 + visualise("SAN_A02_2014", "FirstDegradation", 30, first_degradation.where(first_degradation > 1))
).cols(1)

#get_jrc_data(*sample[["lon_min", "lat_min", "lon_max", "lat_max"]]).rio.reproject(WEBMERCATOR)

Project biomass and jrc to same raster

#deforestations = [deforestation.rio.reproject_match(chm_mean) for chm_mean in chm_means]
#degradations = [degradation.rio.reproject_match(chm_mean) for chm_mean in chm_means]

Mask pixels which have been deforested X time prior

(do I need to account for disturbance as well? Check JRC user guide)

#years_since_degradation = [(sample.year - degradation.where(degradation > 0)) for sample, deforestation in zip(samples, deforestations)]
#years_since_deforestation = [(sample.year - deforestation.where(deforestation > 0)) for sample, deforestation in zip(samples, deforestations)]
#years_since_deforestation = [year.where(year > 0) for year in years_since_deforestation]

chm_means[0].rio.rep

CRS.from_epsg(32721)

import xarray as xr
from tqdm.autonotebook import tqdm

<ipython-input-46-1d2b4d910b30>:2: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  from tqdm.autonotebook import tqdm

survey_name = "SAN_A02_2014"
survey = paisagens_master.loc[survey_name]
desired_resolution = 50
relevant_datasets = ("TransitionMap_MainClasses", "DeforestationYear", "DegradationYear",
                    "LastDeforestationYear", "LastDegradationYear")

# Load canopy height model (CHM)
chm = load_raster(survey.name, "canopy_height_model", 1, from_crs=survey.crs, to_crs=None)

# Aggregate CHM to desired resolution
chm_mean = chm.rio.reproject(chm.rio.crs, resolution=desired_resolution, resampling=Resampling.average)

# Extract relevant JRC datasets
jrc_data = {dataset: get_jrc_data(survey.lon_min, 
                                  survey.lat_min, 
                                  survey.lon_max, 
                                  survey.lat_max, 
                                  dataset=dataset
                                 ).squeeze()#.rio.reproject(survey.crs)
           for dataset in relevant_datasets}

annual_change = xr.concat([get_jrc_data(survey.lon_min, 
                                        survey.lat_min, 
                                        survey.lon_max, 
                                        survey.lat_max, 
                                        dataset="AnnualChange",
                                        year = year
                                       ).squeeze().assign_coords({"year": year})
                             for year in tqdm(range(1992, survey.survey_year + 1))], 
                          dim="year")

# Get bounds
l,b,r,t = gpd.GeoSeries(survey.geometry_bbox, crs=SIRGAS_BRAZIL).to_crs(WEBMERCATOR).bounds.values[0]

xs, ys = np.array(list(gpd.GeoSeries(survey.geometry_metrics, crs=SIRGAS_BRAZIL).to_crs(WEBMERCATOR)[0].exterior.coords)).T
poly = hv.Polygons({"x": xs, "y": ys}).redim(x="Easting", y="Northing").opts(alpha=0.5)

from src.visualisation.visualisation_params import *
from bokeh.models import FixedTicker
from src.processing.jrc_processing import last_of_value

def last_observation(annual_change, jrc_class_value, first_observation=None):
    
    if first_observation is not None:
        last_observation = xr.zeros_like(first_observation)
    else:
        last_observation = xr.zeros_like(annual_change[0])

    # Compute last year of observation (Note: Due to the size of the datasets this takes
    #  ~10 min, even with numba speedup)
    last_observation.data = last_of_value(
        annual_change.data, value=jrc_class_value
    )

    # Fix all negative/fill values and all values prior to 1990 with the first year of
    #  deforestation.
    if first_observation is not None:
        last_observation = xr.concat(
            [first_observation, last_observation + 1990], dim="observations"
        ).max(dim="observations")
    else: 
        last_observation = xr.concat(
            [xr.zeros_like(last_observation), last_observation + 1990], dim="observations"
        ).max(dim="observations")
    
    return last_observation.astype(float)

last_deforested = last_observation(annual_change, 3).rio.reproject(WEBMERCATOR)

last_deforested

<xarray.DataArray (y: 122, x: 232)>
array([[   0.,    0.,    0., ..., 2012., 2012., 2000.],
       [   0.,    0.,    0., ..., 2012., 2012., 2012.],
       [   0.,    0.,    0., ..., 2012., 2012., 2012.],
       ...,
       [2012., 2012., 2012., ...,    0., 2012., 2012.],
       [2012., 2012., 2012., ..., 2012., 2012., 2012.],
       [2012., 2012., 2012., ..., 2012., 2012., 2012.]])
Coordinates:
  * x            (x) float64 -6.113e+06 -6.113e+06 ... -6.106e+06 -6.106e+06
  * y            (y) float64 -3.652e+05 -3.652e+05 ... -3.688e+05 -3.688e+05
    year         int64 1992
    band         int64 1
    spatial_ref  int64 0
Attributes:
    _FillValue:    -9999.0
    grid_mapping:  spatial_ref

mpl.cm.get_cmap('viridis', 30)

<matplotlib.colors.ListedColormap at 0x7ff00938e340>

viridisBig = mpl.cm.get_cmap('viridis', 512)
newcmp = ListedColormap(viridisBig(np.linspace(0.25, 0.75, 256)))
plot_examples([viridis, newcmp])

<matplotlib.colors.ListedColormap at 0x7ff01019c070>

(
hv.Image(annual_change.loc[survey.survey_year].rio.reproject(WEBMERCATOR)).opts(tools=["hover"], 
                                                     **opt_annual_change,
                                                     **opt_size) * poly
+ 
hv.Image(annual_change.where(annual_change==4).loc[survey.survey_year].rio.reproject(WEBMERCATOR)).opts(tools=["hover"], 
                                                     **opt_annual_change,
                                                     **opt_size) * poly
+ 
hv.Image(last_deforested.where(last_deforested > 0), vdims=["deforestation_year"]).opts(**opt_size, 
                                                          colorbar=True,
                                                          cmap= mpl.cm.get_cmap('viridis', 30),
                                                          tools=["hover"]) * poly
).cols(1)

FixedTicker()

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-179-82ebe2fe7f42> in <module>
----> 1 FixedTicker({key: val["name"] for key, val in annual_change_legend.items()})

TypeError: __init__() takes 1 positional argument but 2 were given

def pixel_based_recovery(survey_name)

last_deforestation.rio.reproject_match()

(30.005401194024046, -30.005401194024046)

last_deforestation.rio.clip_box()

was_deforested = last_deforestation > 0
is_recovering

%%time

years_since_deforestation = []
for observation, chm_mean in zip(samples, chm_means):
    print(observation.survey_year)
    # Selectors
    deforested = deforestation > 0
    never_degraded = degradation == 0
    not_distrubed_since_deforested = deforestation > degradation
    deforested_before_observation = deforestation < observation.survey_year
    
    # Convert deforestation year into recovery period
    recovery_period = (observation.survey_year 
                       - deforestation.where(deforested 
                                             & deforested_before_observation 
                                             & (never_degraded | not_distrubed_since_deforested)
                                            )
                      )
    
    # Reproject and match recovery period area to area with lidar observation
    recovery_period = recovery_period.rio.reproject_match(chm_mean)
    years_since_deforestation.append(recovery_period)

2014
2017
2016
CPU times: user 421 ms, sys: 20 ms, total: 441 ms
Wall time: 447 ms

accumulated_heights = [chm_mean.where(recovery_period > 0) for chm_mean, recovery_period in zip(chm_means, years_since_deforestation)]
accumulated_biomasses_longo2016 = [biomass.where(recovery_period > 0) for biomass, recovery_period in zip(longo2016_acds, years_since_deforestation)]
accumulated_biomasses_asner2014 = [biomass.where(recovery_period > 0) for biomass, recovery_period in zip(asner2014_acds, years_since_deforestation)]
accumulated_biomasses_asner2014c = [biomass.where(recovery_period > 0) for biomass, recovery_period in zip(asner2014_colombia_acds, years_since_deforestation)]

xs = np.concatenate([recovery_period.squeeze().data.flatten() for recovery_period in years_since_deforestation])
ys_height = np.concatenate([height.squeeze().data.flatten() for height in accumulated_heights])
ys_longo2016 = np.concatenate([biomass.squeeze().data.flatten() for biomass in accumulated_biomasses_longo2016])
ys_asner2014 = np.concatenate([biomass.squeeze().data.flatten() for biomass in accumulated_biomasses_asner2014])
ys_asner2014c = np.concatenate([biomass.squeeze().data.flatten() for biomass in accumulated_biomasses_asner2014c])

# Visualisation imports
from holoviews import opts
import holoviews as hv
from holoviews.operation.datashader import rasterize
from holoviews.element import tiles
hv.extension("bokeh")

df = pd.DataFrame(np.vstack((xs, 
                        ys_height, 
                        ys_longo2016, 
                        ys_asner2014, 
                        ys_asner2014c)
                      ).T, 
             columns=["recovery_period [years]", 
                      "mean_height [m]", 
                      "longo2016_biomass [kgC/m2]", 
                      "asner2014_biomass [kgC/m2]", 
                      "asner2014c_biomass [kgC/m2]"]
            ).dropna().reset_index(drop=True)
table = hv.Table(df)
points = hv.Points(df)
df.head()

	recovery_period [years]	mean_height [m]	longo2016_biomass [kgC/m2]	asner2014_biomass [kgC/m2]	asner2014c_biomass [kgC/m2]
0	12.0	15.705148	6.877455	9.425895	6.906647
1	12.0	14.114010	5.698736	8.514473	6.031612
2	17.0	17.577612	8.385388	10.492824	7.967184
3	17.0	14.914895	6.280082	8.973813	6.468921
4	12.0	11.979083	4.269915	7.283663	4.898956

hv.BoxWhisker(table, kdims=["recovery_period [years]"], vdims="mean_height [m]").sort().opts(**opt_size)

BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('index', 51), ('mean_height [m]', 51), ('mean_height_left_square_bracket_m_right_square_bracket', 0)

max_observation = int(points["recovery_period [years]"].max())

hist = points.hist(dimension=["recovery_period [years]"], 
            num_bins=max_observation, 
            bin_range=(0.5, max_observation + 0.5),
            adjoin=False).opts(**opt_size, tools=["hover"])
hist

hv.BoxWhisker(table, kdims=["recovery_period [years]"], vdims="mean_height [m]").sort().opts(**opt_size)

BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('index', 20), ('mean_height [m]', 20), ('mean_height_left_square_bracket_m_right_square_bracket', 0)

points = hv.Points(df)

max_observation = int(points["recovery_period [years]"].max())

hist = points.hist(dimension=["recovery_period [years]"], 
            num_bins=max_observation, 
            bin_range=(0.5, max_observation + 0.5),
            adjoin=False).opts(**opt_size, tools=["hover"])
hist

hv.BoxWhisker(table, kdims=["recovery_period [years]"], vdims="longo2016_biomass [kgC/m2]").sort().opts(**opt_size)

BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('index', 58), ('longo2016_biomass [kgC/m2]', 58), ('longo2016_biomass_left_square_bracket_kgC_over_m2_right_square_bracket', 0)

rasterize(points).opts(cmap="kbc_r", cnorm="eq_hist").opts(**opt_size,  tools=["hover"])

Correcting for "last year of deforestation" rather than "first year"

In the JRC TMF data specification it says that the DeforestationYear marks the first year of deforestation. What we need for our regrowth analysis, however, is the last year of deforestation. This will tell us from when on the forest was undistrubed and could resume its growth.

JRC procedure: The JRC-TMF classification process starts out by mapping disturbances in the forest canopy on a yearly basis, regardless of their permanence. The distinction between deforestation and forest degradation is made three years after the disturbance occurred by measuring the permanence of the forest disturbance over time. If the forest canopy is disturbed permanently, i.e. shows no signs of forest regrowth over the three years following the disturbance, the ‘forest disturbance’ pixel falls into the deforestation class. If a ‘forest disturbance’ pixel shows clear signs of forest regrowth within the three years following the disturbance, it is classified as forest degradation. In consequence, the distribution of yearly deforestation and forest degradation areas within the measured yearly overall forest disturbance areas are consolidated until 2017, but are estimated (indicated by stars in Figure 21) for the years 2018-2020 on basis of the 16-year average for the period 2002-2017.

The dataset contains yearly information about disturbances in tropical humid forest, starting from year 1984, and about the trajectories over time of each disturbed forest pixel. It differentiates between short-duration and long-duration disturbances, depending on the forest regrowth after the disturbance. Forest regrowth classes only apply to previously deforested areas, not for forest degradation areas. The dataset differentiates between old, young and very young regrowth, depending on when the forest started to grow back after deforestation. Deforestation and forest degradation information are aggregated to the periods ‘before 2000’, ‘2000-2009’ and ‘2010- 2019’ in order to make the map easier to read. However, the data contains the full information, i.e. the precise years of first or consecutive forest disturbances, plus their classification into the deforestation or forest degradation classes

Note: Text from here

from src.processing.jrc_processing import last_of_value, first_of_value

annual_change_maps = xr.concat(objs=[get_jrc_data(lon_min, 
                                             lat_min, 
                                             lon_max, 
                                             lat_max, 
                                             dataset="AnnualChange", 
                                             year=year
                                            ).squeeze().assign_coords({"year": year}) 
                                      for year in range(1990, 2020)],
                                dim="year")

%%time
deforestation = rxr.open_rasterio("/gws/nopw/j04/forecol/data/JRC/DeforestationYear/JRC_TMF_DeforestationYear_v1_1982_2019_SAM_ID80_N30_W80.tif").squeeze()

CPU times: user 102 ms, sys: 18.7 ms, total: 121 ms
Wall time: 123 ms

annual_change_maps = xr.concat(objs=[rxr.open_rasterio(
                        f"/gws/nopw/j04/forecol/data/JRC/AnnualChange/tifs/JRC_TMF_AnnualChange_v1_{year}_SAM_ID80_N30_W80.tif"
                                    ).squeeze().assign_coords({"year": year}) 
                                      for year in range(1990, 2020)],
                                dim="year")

%%time
last_deforestation = xr.zeros_like(deforestation)

CPU times: user 23.2 s, sys: 17 s, total: 40.2 s
Wall time: 41.5 s

%%time
last_deforestation.data = last_of_value(annual_change_maps.data, 3)

CPU times: user 9min 51s, sys: 47.1 s, total: 10min 38s
Wall time: 10min 38s

%%time
last_deforestation = xr.concat([deforestation.squeeze(), last_deforestation + 1990], dim="observations").max(dim="observations")

CPU times: user 46.2 s, sys: 53.3 s, total: 1min 39s
Wall time: 1min 45s

from src.constants import JRC_PATH

%%time
last_deforestation.rio.to_raster(JRC_PATH / "LastDeforestationYear" / "JRC_TMF_LastDeforestationYear_v1_1982_2019_SAM_ID80_N30_W80.tif",
                                compress='lzw')

CPU times: user 56.1 s, sys: 28.6 s, total: 1min 24s
Wall time: 1min 31s

%%time
last_deforestation = xr.zeros_like(deforestation.squeeze())
last_deforestation.data = last_of_value(annual_change_maps.data, 3)
last_deforestation = xr.concat([deforestation.squeeze(), last_deforestation + 1990], dim="observations").max(dim="observations")

<xarray.DataArray (year: 30, y: 37459, x: 37108)>
array([[[5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        ...,
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5]],

       [[5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        ...,
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5]],

       [[5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        ...,
...
        ...,
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5]],

       [[5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        ...,
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5]],

       [[5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        ...,
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5]]], dtype=uint8)
Coordinates:
    band         int64 1
  * y            (y) float64 30.09 30.09 30.09 30.09 ... 20.0 20.0 20.0 20.0
  * x            (x) float64 -80.0 -80.0 -80.0 -80.0 ... -70.0 -70.0 -70.0 -70.0
    spatial_ref  int64 0
  * year         (year) int64 1990 1991 1992 1993 1994 ... 2016 2017 2018 2019
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    long_name:     Dec1990
    grid_mapping:  spatial_ref

rxr.open_rasterio("/gws/nopw/j04/forecol/data/JRC/AnnualChange/tifs/JRC_TMF_AnnualChange_v1_2019_SAM_ID80_N30_W80.tif")

<xarray.DataArray (band: 1, y: 37459, x: 37108)>
[1390028572 values with dtype=uint8]
Coordinates:
  * band         (band) int64 1
  * y            (y) float64 30.09 30.09 30.09 30.09 ... 20.0 20.0 20.0 20.0
  * x            (x) float64 -80.0 -80.0 -80.0 -80.0 ... -70.0 -70.0 -70.0 -70.0
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    long_name:     Dec2019
    grid_mapping:  spatial_ref

i = 0
y, x = regrowth_pixels[i]
regrowth_pixels = np.argwhere(transition_map.data.squeeze() == 30)  # 30 = regrowth

# Now I need to get the last deforestation year
first_deforested = np.argwhere(annual_change_maps.data[:, y, x] == 3).min()
print(f"First deforested: {1990 + first_deforested}")

# Now I need to get the last deforestation year
last_deforested = np.argwhere(annual_change_maps.data[:, y, x] == 3).max()
print(f"Last deforested: {1990 + last_deforested}")

First deforested: 1995
Last deforested: 2014

%%time
last_deforestation = xr.zeros_like(deforestation.squeeze())
last_deforestation.data = last_of_value(annual_change_maps.data, 3)
last_deforestation = xr.concat([deforestation.squeeze(), last_deforestation + 1990], dim="observations").max(dim="observations")

CPU times: user 39.1 ms, sys: 61 µs, total: 39.2 ms
Wall time: 38.7 ms

np.sum(deforestation > 0)

<xarray.DataArray ()>
array(19750)
Coordinates:
    band         int64 1
    spatial_ref  int64 0

np.sum(last_deforestation > 0)

<xarray.DataArray ()>
array(19750)
Coordinates:
    band         int64 1
    spatial_ref  int64 0

%%time
first_of_value(annual_change_maps.data, 3)

CPU times: user 32.7 ms, sys: 4.63 ms, total: 37.3 ms
Wall time: 37.3 ms

array([[-9999,    26, -9999, ...,     8,     8,    14],
       [-9999,    26, -9999, ...,     8,     8,    14],
       [-9999,    26, -9999, ...,    14,    14,    14],
       ...,
       [-9999, -9999, -9999, ..., -9999, -9999, -9999],
       [-9999, -9999, -9999, ..., -9999, -9999, -9999],
       [-9999, -9999, -9999, ..., -9999, -9999, -9999]])

%%time
last_deforested_year = xr.zeros_like(transition_map.squeeze(), dtype=np.uint16)
for y in range(len(annual_change_maps.y)):
    for x in range(len(annual_change_maps.x)):
        deforested_years = np.argwhere(annual_change_maps.data[:, y, x] == 3)
        if len(deforested_years) > 0:
            last_deforested_year[y, x] = 1990 + deforested_years.max()
            
last_deforested_year = xr.concat([last_deforested_year, deforestation.squeeze()], dim="disturbance").max(axis=0)

CPU times: user 3.77 s, sys: 141 ms, total: 3.91 s
Wall time: 3.85 s

np.argwhere(annual_change_maps.where(annual_change_maps == 3)[:, y:y+3 , x:x+3].data == 3)

array([[ 5,  0,  0],
       [ 5,  0,  1],
       [ 5,  0,  2],
       [ 5,  1,  2],
       [ 5,  2,  0],
       [ 5,  2,  1],
       [ 5,  2,  2],
       [ 6,  0,  0],
       [ 6,  0,  1],
       [ 6,  0,  2],
       [ 6,  1,  2],
       [ 6,  2,  0],
       [ 6,  2,  1],
       [ 6,  2,  2],
       [ 7,  0,  0],
       [ 7,  0,  1],
       [ 7,  0,  2],
       [ 7,  1,  2],
       [ 7,  2,  0],
       [ 7,  2,  1],
       [ 7,  2,  2],
       [ 8,  0,  0],
       [ 8,  0,  1],
       [ 8,  0,  2],
       [ 8,  1,  2],
       [ 8,  2,  0],
       [ 8,  2,  1],
       [ 8,  2,  2],
       [ 9,  0,  0],
       [ 9,  0,  1],
       [ 9,  0,  2],
       [ 9,  1,  2],
       [ 9,  2,  0],
       [ 9,  2,  1],
       [ 9,  2,  2],
       [10,  0,  0],
       [10,  0,  1],
       [10,  0,  2],
       [10,  1,  2],
       [10,  2,  0],
       [10,  2,  1],
       [10,  2,  2],
       [11,  0,  0],
       [11,  0,  1],
       [11,  0,  2],
       [11,  1,  2],
       [11,  2,  0],
       [11,  2,  1],
       [11,  2,  2],
       [12,  0,  0],
       [12,  0,  1],
       [12,  0,  2],
       [12,  1,  2],
       [12,  2,  0],
       [12,  2,  1],
       [12,  2,  2],
       [13,  0,  0],
       [13,  0,  1],
       [13,  0,  2],
       [13,  1,  2],
       [13,  2,  0],
       [13,  2,  1],
       [13,  2,  2],
       [14,  0,  0],
       [14,  0,  1],
       [14,  0,  2],
       [14,  1,  2],
       [14,  2,  0],
       [14,  2,  1],
       [14,  2,  2],
       [15,  0,  0],
       [15,  0,  1],
       [15,  0,  2],
       [15,  1,  2],
       [15,  2,  0],
       [15,  2,  1],
       [15,  2,  2],
       [16,  0,  0],
       [16,  0,  1],
       [16,  0,  2],
       [16,  1,  2],
       [16,  2,  0],
       [16,  2,  1],
       [16,  2,  2],
       [17,  0,  0],
       [17,  0,  1],
       [17,  0,  2],
       [17,  1,  2],
       [17,  2,  0],
       [17,  2,  1],
       [17,  2,  2],
       [18,  0,  0],
       [18,  0,  1],
       [18,  0,  2],
       [18,  1,  2],
       [18,  2,  0],
       [18,  2,  1],
       [18,  2,  2],
       [19,  0,  0],
       [19,  0,  1],
       [19,  0,  2],
       [19,  1,  2],
       [19,  2,  0],
       [19,  2,  1],
       [19,  2,  2],
       [20,  0,  0],
       [20,  0,  1],
       [20,  0,  2],
       [20,  1,  2],
       [20,  2,  0],
       [20,  2,  1],
       [20,  2,  2],
       [21,  0,  0],
       [21,  0,  1],
       [21,  0,  2],
       [21,  1,  2],
       [21,  2,  0],
       [21,  2,  1],
       [21,  2,  2],
       [22,  0,  0],
       [22,  0,  1],
       [22,  0,  2],
       [22,  1,  2],
       [22,  2,  0],
       [22,  2,  1],
       [22,  2,  2],
       [23,  0,  0],
       [23,  0,  1],
       [23,  0,  2],
       [23,  1,  2],
       [23,  2,  0],
       [23,  2,  1],
       [23,  2,  2],
       [24,  0,  0],
       [24,  0,  1],
       [24,  0,  2],
       [24,  1,  2],
       [24,  2,  0],
       [24,  2,  1],
       [24,  2,  2],
       [25,  0,  1],
       [25,  0,  2],
       [25,  1,  2],
       [25,  2,  0],
       [25,  2,  1],
       [25,  2,  2],
       [26,  0,  1],
       [26,  0,  2],
       [26,  1,  2],
       [26,  2,  0],
       [26,  2,  1],
       [26,  2,  2],
       [27,  0,  1],
       [27,  0,  2],
       [27,  1,  2],
       [27,  2,  0],
       [27,  2,  1],
       [27,  2,  2],
       [28,  0,  1],
       [28,  0,  2],
       [28,  1,  2],
       [28,  2,  0],
       [28,  2,  1],
       [28,  2,  2]])

np.argwhere(annual_change_maps.where(annual_change_maps == 3)[:, y:y+10 , x:x+10].data == 3).max(axis=1).shape

(2331,)

import numba

deforestation.squeeze()

<xarray.DataArray (y: 161, x: 462)>
array([[   0., 2016.,    0., ..., 1998., 1998., 2004.],
       [   0., 2016.,    0., ..., 1998., 1998., 2004.],
       [   0., 2016.,    0., ..., 2004., 2004., 2004.],
       ...,
       [   0.,    0.,    0., ...,    0.,    0.,    0.],
       [   0.,    0.,    0., ...,    0.,    0.,    0.],
       [   0.,    0.,    0., ...,    0.,    0.,    0.]])
Coordinates:
    band         int64 1
  * y            (y) float64 -3.279 -3.279 -3.279 -3.28 ... -3.321 -3.322 -3.322
  * x            (x) float64 -54.94 -54.94 -54.94 ... -54.82 -54.82 -54.82
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    grid_mapping:  spatial_ref

last_deforested_year.where(last_deforested_year > 0)

<xarray.DataArray ()>
array(19750)
Coordinates:
    band         int64 1
    spatial_ref  int64 0

(last_deforested_year - deforestation.squeeze()).min()

<xarray.DataArray ()>
array(-1989.)
Coordinates:
    band         int64 1
    spatial_ref  int64 0

np.sum(deforestation > 0)

<xarray.DataArray ()>
array(19750)
Coordinates:
    spatial_ref  int64 0

np.sum(last_deforested_year > 0)

<xarray.DataArray ()>
array(19717)
Coordinates:
    band         int64 1
    spatial_ref  int64 0

<xarray.DataArray (year: 29)>
array([1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
       3, 3, 3, 4, 4, 4, 4], dtype=uint8)
Coordinates:
    band         int64 1
    y            float64 -3.279
    x            float64 -54.93
    spatial_ref  int64 0
  * year         (year) int64 1990 1991 1992 1993 1994 ... 2015 2016 2017 2018
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    long_name:     Dec1990
    grid_mapping:  spatial_ref

deforestation[:, y,x]

<xarray.DataArray (band: 1)>
array([1995.])
Coordinates:
  * band         (band) int64 1
    y            float64 -3.279
    x            float64 -54.93
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    grid_mapping:  spatial_ref

degradation[:, y,x]

<xarray.DataArray (band: 1)>
array([0.])
Coordinates:
  * band         (band) int64 1
    y            float64 -3.279
    x            float64 -54.93
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    grid_mapping:  spatial_ref

annual_change_maps[:, 1, -2]

<xarray.DataArray (year: 29)>
array([1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
       3, 3, 3, 3, 3, 3, 3], dtype=uint8)
Coordinates:
    band         int64 1
    y            float64 -3.279
    x            float64 -54.82
    spatial_ref  int64 0
  * year         (year) int64 1990 1991 1992 1993 1994 ... 2015 2016 2017 2018
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    long_name:     Dec1990
    grid_mapping:  spatial_ref

transition_map[:, 2,1] # 30 = forest regrowth

<xarray.DataArray (band: 1)>
array([30], dtype=uint8)
Coordinates:
  * band         (band) int64 1
    y            float64 -3.279
    x            float64 -54.94
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    long_name:     b1
    grid_mapping:  spatial_ref

deforestation[:, 2, 1]

<xarray.DataArray (band: 1)>
array([2016.])
Coordinates:
  * band         (band) int64 1
    y            float64 -3.279
    x            float64 -54.94
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    grid_mapping:  spatial_ref

degradation[:,2,1]

<xarray.DataArray (band: 1)>
array([2008.])
Coordinates:
  * band         (band) int64 1
    y            float64 -3.279
    x            float64 -54.94
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    grid_mapping:  spatial_ref

annual_change_maps[:,2, 1]

<xarray.DataArray (year: 29)>
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2,
       2, 2, 2, 2, 3, 4, 4], dtype=uint8)
Coordinates:
    band         int64 1
    y            float64 -3.279
    x            float64 -54.94
    spatial_ref  int64 0
  * year         (year) int64 1990 1991 1992 1993 1994 ... 2015 2016 2017 2018
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    long_name:     Dec1990
    grid_mapping:  spatial_ref

annual_change_maps[:,0,0]

<xarray.DataArray (year: 28)>
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 2, 2], dtype=uint8)
Coordinates:
    band         int64 1
    y            float64 -3.279
    x            float64 -54.94
    spatial_ref  int64 0
  * year         (year) int64 1990 1991 1992 1993 1994 ... 2014 2015 2016 2017
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    long_name:     Dec1990
    grid_mapping:  spatial_ref

Legacy

Extract survey names

#survey_year = pd.Series(paisagens_master.index, index=paisagens_master.index).str.findall("\w+_(\d{4})\w*").apply(lambda x: int(x[0]))
#paisagens_master["survey_year"] = survey_year
#eba_master["survey_year"] = eba_master.year

#paisagens_master.drop(columns=["geometry_metrics", "geometry_bbox"]).to_csv(PAISAGENSLIDAR_PATH / "paisagens_master_table.csv")
#gpd.GeoDataFrame(paisagens_master, crs=SIRGAS_BRAZIL, geometry="geometry_metrics").to_feather(PAISAGENSLIDAR_PATH / "paisagens_master_table.feather")
#gpd.GeoDataFrame(paisagens_master, crs=SIRGAS_BRAZIL, geometry="geometry_metrics").to_parquet(PAISAGENSLIDAR_PATH / "paisagens_master_table.parquet")

#eba_master.drop(columns=["geometry_metrics", "geometry_bbox", "geometry_metadata"]).to_csv(EBALIDAR_PATH / "eba_master_table.csv")
#gpd.GeoDataFrame(eba_master, crs=SIRGAS_BRAZIL, geometry="geometry_metrics").to_feather(EBALIDAR_PATH / "eba_master_table.feather")
#gpd.GeoDataFrame(eba_master, crs=SIRGAS_BRAZIL, geometry="geometry_metrics").to_parquet(EBALIDAR_PATH / "eba_master_table.parquet")

Extract latitude longitude

e_lon_lat = gpd.GeoSeries(eba_master.geometry_bbox, 
              crs=SIRGAS_BRAZIL).to_crs(WGS84).bounds.rename(columns={"minx": "min_lon",
                                                               "miny": "min_lat", 
                                                               "maxx": "max_lon",
                                                               "maxy": "max_lat"})

p_lon_lat = gpd.GeoSeries(paisagens_master.geometry_bbox, 
              crs=SIRGAS_BRAZIL).to_crs(WGS84).bounds.rename(columns={"minx": "min_lon",
                                                               "miny": "min_lat", 
                                                               "maxx": "max_lon",
                                                               "maxy": "max_lat"})